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We investigate the optimal control of open quantum systems, in particular, the mutual influence 
of driving and dissipation. A stochastic approach to open-system control is developed, using a 
generalized version of Krotov's iterative algorithm, with no need for Markovian or rotating-wave 
approximations. The application to a harmonic degree of freedom reveals cooperative effects of 
driving and dissipation that a standard Markovian treatment cannot capture. Remarkably, control 
can modify the open-system dynamics to the point where the entropy change turns negative, thus 
achieving cooling of translational motion without any reliance on internal degrees of freedom. 
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Introduction. The control of quantum dynamics for the 
accurate preparation of a prescribed quantum state by a 
tailored time-dependent field is a task of key importance 
in quantum physics and related disciplines. With increas- 
ing complexity of devices for quantum information pro- 
cessing the destructive role of environmental fluctuations 
has become a severe limitation to further progress. For 
example, neutral atoms or ions in electromagnetic traps 
are exposed to fluctuations of (comparatively hot) chip 
surfaces [I], while in superconducting circuits diffusing 
charges and electromagnetic fluctuations affect fidelities 
quite substantially [2J. 

In this context, optimal control theory (OCT) has 
emerged as a key ingredient in strategies to tame the 
effect of decoherence and other imperfections either di- 
rectly by mitigating their effect or indirectly by speeding 
up operations. In limiting cases or under idealized condi- 
tions OCT has already been used to increase the fidelity 
of simple quantum gates [3H7] or to construct more com- 
plex protocols [S]. Generally, however, OCT has treated 
environmental interactions mostly by heuristic or approx- 
imate methods so far. A simple strategy has been fol- 
lowed in Ref. [S], where the impact of a heat bath was 
taken into account by assuming an initial thermal state 
while neglecting its effect completely during its dynamics. 
Some progress has been achieved within fully dynamical 
approaches based on standard Markovian master equa- 
tions. Unfortunately, these schemes have severe draw- 
backs: First, they become inconsistent for strong control 
fields unless additional field-dependent memory terms in 
the dissipator are introduced [TO] HI] . Second, the ro- 
tating wave approximation (RWA), which is usually per- 
formed on the system-environment interaction in order 
to obtain a master equation with complete positivity, is 
known to be unreliable in driven systems |10j . Third, at 
sufficiently low temperatures or for reservoirs with struc- 
tured spectral mode densities the true dynamics of open 
quantum systems is non-Markovian. These errors in the 
dynamics may be further amplified by OCT search algo- 



rithms. As a result, OCT computations using standard 
master equations yield mismatched fields, which perform 
poorly when applied to a realistic setting. 

In this Letter we present a treatment of optimal con- 
trol of open quantum systems which is not susceptible 
to these problems. Stochastic Liouville-von Neumann 
(SLN) equations (T2J02] provide an approach to quantum 
dissipation in a driven system which is both conceptually 
transparent and formally exact. Field-induced modifica- 
tions of environmental effects are thus included a priori. 
From the SLN equations one arrives naturally at the def- 
inition of a state-costate pair [13] of dynamical variables 
with the simple first-order equations of motion required 
by gradient-based and related OCT methods. 

We outline the salient features of our OCT technique 
(for details see Supplemental Material HS]) and apply 
it to a simple model. Notably, it turns out that con- 
trol fields extracted from a RWA-based master equation 
differ substantially from exact ones obtained within the 
SLN scheme. This leads even to the counterintuitive phe- 
nomenon of control-induced cooling, which is completely 
missing in the RWA approach. 

Open-system dynamics and control algorithm. — We 
start with the exact stochastic equation of motion |12j 

+ ~Z{t)[q,Qt, v {t)] + ±v(t){q,Qt, v } (1) 

for noisy samples Q£, v {t) of the density matrix of an 
open quantum system. The system is characterized by a 
Hamiltonian H s governing its autonomous motion, aug- 
mented by a term H c (t) describing the influence of time- 
dependent control fields Uj(t). It interacts bilinear ly 
with a dissipative environment whose quantum statis- 
tical force-force correlation function is mapped to the 
noise statistics as follows [12] : (i) the autocorrelation of £ 
matches the quantum noise of the reservoir, (ii) the cross- 
correlations of £ and v match the dynamical response 
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of the environment and (iii) the correlations (v(t)v{t')) 
vanish identically. According to the last condition v{t) 
is a complex variable with random phase. The non- 
vanishing correlations are identical to the real and imag- 
inary parts of the kernel defining the Feynman- Vernon 
influence functional [Jj)J E] ; from which Eq. ([!]) can be 
derived without approximations. 

The fluctuations and the dynamical response of the 
thermal environment are fully characterized by its tem- 
perature and spectral density J(uj), which in the present 
case will be taken as Ohmic (proportional to uj up to a 
UV cutoff u c ). Both 1/uj c and, more important, the ther- 
mal time H/3 are intrinsic time scales of the environmental 
fluctuations. Eq. now simplifies to 

~^o[q, {p, «(*)}] + «(*)] . (2) 

where 70 is the damping rate of a Brownian particle of 
mass m |18j . The physical density matrix is a stochastic 
average of the form g(t) = E[g^(t)]. 

At the price of introducing an explicit noise variable 
Eq. ([2]) represents the exact non-Markovian dy- 
namics in terms of a stochastic ensemble with time-local 
equations of motion. All memory effects inherent in the 
reservoir dynamics are contained in the time-dependent 
correlations. There is no decomposition of the envi- 
ronment into additional degrees of freedom and a sec- 
ondary, Markovian environment [19) . In an extreme high- 
temperature limit, Eq. |2]) becomes Markovian and re- 
duces to the master equation of Caldeira and Leggett [20] . 

Optimal control means searching for control signals 
which drive desirable characteristics of the dynamics to 
extremal values. Here we consider optimization objec- 
tives defined by minimization of an expectation value 
(M) at a specified end time T. This leads to a search 
for extrema of the objective functional 

B[u(i)] = E[tr{Mg£(T)}] = tv{Mg(T)}, (3) 

where the first equality relies on the map u i-> g$(T) 
implicit in the initial- value problem of Eq. ^ . Alterna- 
tively, Eq. Q can be interpreted as a constraint on si- 
multaneous variations of u(t) and g^(t). This constraint 
needs to be taken into account through a time-dependent 
Lagrange multiplier A^(t), which also depends on the par- 
ticular noise realization Variational calculus leads 
to the equation of motion 

A*(t) = -£%(t) = ~[H B) At(t)] ~ £[ffc,M*)] 

-± JO { P)[qi A^t)}} + ^(t)[Q,M(t)} (4) 

for A^(t), called the costate of the optimal control prob- 
lem in this context. Equation Q is not an initial value 



problem; it needs to be solved with the terminal bound- 
ary condition A^(T) + M = arising from variation of 
the final state. 

Now the gradient of the objective functional under the 
above constraint is given by 
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where g^(t) and Aj(i) obey the stochastic equations of 
motion (2j) and Q. 

The preceding considerations show that the SLN ap- 
proach treats quantum memory effects in a mathemati- 
cal language which integrates seamlessly into the state- 
cost ate framework of standard OCT techniques. For nu- 
merical computations, we have adapted the monoton- 
ically convergent algorithm of Krotov [3TJ [52] to the 
present case of non-Markovian stochastic propagation. 
This algorithm improves performance by substituting 
gradient search steps with a nonlocal generalization of 
Eq. ([5]). Key performance characteristics are improved 
by this change (see |15j). 

Application. — As a common model we consider a har- 
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monic oscillator, i.e., H s = ^— + ^^-q 2 , which is subject 
to an additional potential H c (t) = —F(t)q + ^/A{t)q 2 
depending on a force control fields F(t) = Ui(t) and a 
tuning control field A(£) = 1*2 (i), which modifies the in- 
stantaneous resonance frequency, u 2 = Uq + A. This 
choice is not only a model for typical realizations as, e.g., 
trapped atoms or ions or low energy dynamics of Joseph- 
son junctions, it also offers itself as a simple test case 
where an intuitive interpretation of numerical results may 
be feasible. It is nontrivial since it includes the nonlinear 
response of the system to parametric driving, and it is 
fairly generic as it applies to potentials with harmonic 
minima. 

Under the equation of motion ([2]), the individual sam- 
ples g^ remain Gaussian for Gaussian initial states. This 
allows us to rephrase the equation of motion ^ as a sys- 
tem of ordinary stochastic differential equations for the 
first and second cumulants (means and variances) associ- 
ated with g$, i.e., (q) c , (p) c , (q 2 ) c , (p 2 ) c and {\{qp+pq)) c - 
A similar consideration holds for the costate dynamics 
Q if a maximal overlap with a Gaussian target state 
is chosen as optimization objective, i.e., M = 1 — A, 
where A = |a)(a| projects onto a coherent state. We 
thus obtain closed equations of motion for the first two 
cumulants in the propagation of both the state g^{t) and 
the costate A^(t). While the effect of the linear control 
F(t) alone is given by linear response theory, the dy- 
namical squeezing through a time-dependent A(t) leads 
to nontrivial dynamics, as does the combined action of 
both controls. We have explored these effects numeri- 
cally, computing the expectation values in Eq. ^ explic- 
itly through a large number of samples (typically 10 ). 
This has the advantage of being securely based on first 
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principles, without resorting to approximations of the dy- 
namics. 

In the following, we use natural units (K = kfe = 
1, units luq for energies, angular frequencies, or 
rates, 1/y'mwo for lengths, and y/mujQ for momenta). 
We choose a minimal-uncertainty wavepacket centered 
around q = 1 and p = as both initial and target state. 
Values of the temperature and the damping constant are 
chosen in the range typical of superconducting solid-state 
devices [2]. The propagation time T = 20 is roughly 
comparable to the relaxation time in the examples to be 
discussed. 

We compare the results of iteratively determined con- 
trol fields for three types of dynamics inserted for state 
and costate in Eq. |5]): (a) SLN dynamics; (b) the stan- 
dard Markovian master equation of the harmonic oscil- 
lator (23] . with the usual raising and lowering operators 
associated with H s as Lindblad operators; (c) quantum 
dynamics without dissipation. 

Figures [I] and [2] show time-frequency signatures of 
the controls F(t) and A(t) obtained through the win- 
dowed Fourier transform (also short-time Fourier trans- 
form, STFT) using a Gaussian window. Both controls 
show marked differences between the SLN and RWA 
cases. The tendency for more pronounced and more com- 
plex high-frequency features in the SLN case indicates 
the importance of exercising control also on time scales 
of the environmental fluctuations (of order j3) , similar to 
a known strong-field approach to the suppression of de- 
coherence known as 'bang-bang control' [24 . A second 
tendency seen in the SLN results is the application of 
fields spread out over the entire time interval, as com- 
pared to the emphasis on a stronger initial perturbation 
in the cases of RWA dissipation or no dissipation. 

Values of the objective functional achieved with the 
SLN fields for different temperatures and damping con- 
stants are compared in Table [I] Free dynamics (no con- 
trol) would result in values roughly equal to 1/2 for all 
parameters listed. A test of the control fields obtained 
in RWA, inserted in the exact equation of motion, typ- 
ically yields values of the objective functional which are 
up to 100% larger than for controls computed using SLN 
dynamics. The algorithmic property of monotonic con- 
vergence is confirmed by our numerical results. 
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0.0072 


0.0133 



Table I. Results for the minimization of tr{Mp(T)} for vari- 
ous inverse temperatures j3 and different damping constants 
70 in the range typical of mesoscopic quantum circuits or 
condensed-phase chemical reactions. 




Figure 1. (color online) Windowed Fourier transform of the 
optimal control force F(t) obtained using different dynamical 
equations: (a) SLN equation Q, SLN, (b) a simple gener- 
alization of the standard master equation to driven systems, 
RWA, and (c) unitary propagation. Parameters are 70 = 0.05, 
lo c = 50, P = 1. 







RWA 






1 O^ 2 


no dissipation 


B 

I 





0.8 
0.6 
0.4 
0.2 



■0.05 

%04 

0.03 

,0.02 

yo.01 

■0.2 

%A5 

0.1 
0.05 



10 



15 



20 



Figure 2. (color online) Windowed Fourier transform of the 
optimal tuning field A = lj 2 — Uq obtained using dynamical 
equations as in Fig. [I] Different color scales apply to the three 
scenarios. 



Dynamical cooling. — Optimal control for closed sys- 
tems conserves entropy like any unitary time evolution. 
Quantum dissipation invariably creates mixed states in 
the subsystem of interest, i.e., if the initial state is pure 
the entropy of the open system will increase. But can 
optimal control of an open system prevent this or even 
lower the entropy in other cases? To investigate this ques- 
tion, we choose the oscillator ground state as target and 
prepare both system and environment as thermal states 
with equal inverse temperature = 1. In this symmet- 
ric setting, the field F(t) is not needed, since it changes 
the position, but not the shape of the wave packet. We 
therefore consider only the control field A(t) in the fol- 
lowing. The von Neumann entropy of the mixed state 

is given by [25] S(g) = g (v V) c(p 2 ) c - (pq + qp)ll *4j 
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Figure 3. (color online) An open quantum system initially 
equilibrated with its surroundings loses entropy S under an 
optimized control field (solid). In contrast, the standard 
Markovian/RWA master equation leads to increased entropy 
under driving (dashed, see [15] 



with g(x) = (x + §) log (x + |) - (x - §) log (x-\). 
We thus obtain the counterintuitive result that a time- 
dependent control field can modify dissipative dynamics 
to the point where its entropy change turns negative (Fig. 
[3]). We attribute this phenomenon to the cooperative ef- 
fect of driving and dissipation, since neither of the two 
alone can cause this. The subsystem energy of the final 
state decreases below its original thermal value, indicat- 
ing a dynamical cooling effect. In contrast, it can be 
shown (see Supplemental Material [TS]) that commonly 
used RWA methods predict heating above the environ- 
mental temperature for non-zero driving. Consistent cor- 
rections of master equations for finite H c prove to be a 
formidable challenge [TTJ. Moreover, even if H c could be 
used in the construction of the dissipator, the distinction 
between co- and counter-rotating terms would hardly be 
justified. If the control fields change on the timescale of 
the reservoir fluctuations, a 'wobbly frame' rather than 
a rotating frame results. 

In contrast to recent proposals for quantum refrigera- 
tors [27] , which rely on intricate band or level struc- 
tures, we have chosen a model with minimal structure. 
The cooling effect found here seems to be a feature of 
temporal patterns, not of a specifically designed system. 
We also note that no internal degree of freedom is needed 
for the effect to occur. 

Conclusions. — The present SLN approach to optimal 
control enjoys two natural advantages compared to con- 
trol theory based on standard Markovian master equa- 
tions: (i) its noise statistics are by construction inde- 
pendent of the quantum dynamics, i.e., strong external 
driving introduces no need for correction terms, and (ii) 
one arrives at the usual state-costate picture required 
by OCT methods in a straightforward way. Numerical 
control of a harmonic degree of freedom is demonstrated 
with varying parameters and objectives. Most results 
show marked differences compared to the RWA approach, 



where the influence of driving on dissipation is neglected. 
Efficient computations are feasible for environmental cou- 
plings from weak damping up to a quality factor as low as 
Q s» 10. This allows applications to solid-state devices 
such as superconducting circuits with Josephson junc- 
tions and condensed-matter phenomena such as reactive 
dynamics of small molecules in a solvent or on a sur- 
face. Optimal control of a dissipative quantum system 
can extract entropy from a system initially at the same 
temperature as its environment. Dynamical cooling in a 
simple system without special structural features may be 
considered as a likely strategy for mesoscopic quantum 
refrigeration. 
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This document serves as a supplement to our article "Optimal control of open quantum 
systems: cooperative effects of driving and dissipation" [Physical Review Letters (2011)]. It 
describes the numerical algorithm used, a generalization of the Krotov iterative algorithm 
to open quantum systems evolving under a stochastic equation of motion. A comparison 
with the Markovian dynamics of a commonly used quantum master equation shows that the 
latter cannot reproduce the cooperative cooling effect described in the main text. 



I. OBJECTIVE AND CONSTRAINT 



We outline a modification of Krotov's algorithm for the optimization of final-state properties of 
an open quantum system governed by the stochastic differential equation 

&(t) = 2g ( (t) := - l -[H s ,^(t)] - l -[H c {t), g 6 (t)} - ^ l0 [q, {p, g^t)}} + ^(t)[$, &(t)] , (1) 

which must be considered as a dynamical constraint in the context of the optimization problem. 
Here 70 is the damping rate of a Brownian particle in the (ohmic) environment, and is a noise 
variable whose statistics are governed by the quantum correlation function of the environmental 
fluctuations. 

The objective functionals to be minimized is assumed to depend linearly on the terminal state, 
i.e., the functional is the expectation value of an observable M at the final time T, typically a 
projection operator. The physical terminal state g(T) is the stochastic average of the dynamical 
state g^(T), i. e., the objective functional is of the form 

B[u(t)] = tv{Mg(T)} = E[tr{M&(r)}] (2) 

where the vector u(t) denotes control fields which define the time-dependence of H c (t), and the 
first equality implies the equation of motion (1). We use the symbol E for the expectation value 
with respect to the noise only, not for the further average given by tracing over g^(t). 
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II. EXTENDED OBJECTIVE FUNCTIONAL AND OPTIMAL CONTROL EQUATIONS 

Following Refs. [20, 21] (see main text), we use an extended objective functional 

L ? [u(i),^(t);^] := Gt(Q 6 (T)) - f dt R^t, u(t), Q ( (t)) - ^(0, &(0)), (3) 

J o 

containing the auxiliary, real-valued function 4>^(t, gAt)), which depends on the time t, the noise 
realization £(t), and the state sample g^(t). Further auxiliary functions 

G ( (Q((T)) := tr{M^(T)} + ^(T,^(T)), (4) 
Kf(t,u(t),^(t), • ) := tr{ • &g ( (t)}, (5) 

Re(t, u(t), &(t)) := k € (t, u(t), &(t), ^^(*, &(*))J + ^(*, &(*)) (6) 

are chosen such that the expectation value E[Lf] of the extended functional coincides with the 
ordinary objective functional for admissible processes, i. e., pairs (u(i), {&(*)}) of a vector of control 
fields and a stochastic ensemble {g^(t)} with the property that the equation of motion (1) is obeyed 
for the given vector u(t) and for all noise realizations of £(t). In this case the partial derivatives in 
the definition of Rj combine to form a total derivative, effectively canceling the boundary terms of 
contained in Eqs. (3) and (4). 

Different choices of 0j modify the objective functional for non-admissible processes, such as 
the combination of control fields and states from different iterations. For Krotov's algorithm it is 
necessary to choose a function (f)^ with the following property: For any admissible process, is at 
a maximum with respect to variations of g^ with the function u(i) kept fixed, i. e., variations away 
from the submanifold of admissible processes decrease Lj. This means that an appropriate choice 
of (j)^ must conform to the two conditions 

G ( (gf(T)) = max{G^(T))}, (7) 
Rt(t,uW(t),$\t)) = min{R^,u«(t),^))}, (8) 

where the superscript (J) denotes the admissible process resulting from an initial guess or a previous 
iteration. 

After making such a choice of <p^\ we demand that, in addition to being admissible, the 
new process (u^ +1 \t), {f?£ (i)}) is taken from the class of pairs (u(i), {§^{t)}) which fulfill the 
condition 



u(t) = argmaxE[R ? (i,u(t),^(t))] 
u(t) 



(9) 
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at any time t. This condition selects a specific pair (u^ +1 \t),{g^ +1 \t)}) from the set of admissible 
process. Using Eq. (9), the update u^\t) — > uV +1 > (t) can be performed locally in conjunction 
with the propagation of the new state g^ +1 ^ (t) since the dependence of the dynamical state on the 
control fields obeys causality. 

What remains to be specified is a mathematical construction which defines specific computa- 
tional steps using the function </>£. For the linear dynamics of Eq. (1), it turns out that the function 
<j>^(t,g^) appears in the conditions (7) and (8) only in form of its derivative with respect to the 
state, and since a complete derivative can be formed from them, the co-state — dcp^/dg^ can be 
characterized as an additional dynamical variable with the equation of motion 

ksit) = -£%(t) = ~[H B Mt)l ~ ^[Hc(t)Mt)\ ~ ^7o{]3, [q,km + A*(*)]. (10) 

Solving this equation with the end-time boundary condition Aj(T) + M = 0, resulting from (7), 
amounts to a (fully sufficient) construction of (f>^ to first order in variations of g^ around gJ . 

In a variant of the algorithm, described by Sklarz and Tannor [21], the implicit condition (9) is 
substituted by an explicit form 



«? +1> (') =«<?(*)+ ' 



A fc (i) 



(") 



which is equivalent to Eq. (9) if a suitable (j-dependent) cost functional is added to B. It is to be 
noted that the simultaneous appearance of j and j + 1 on the r.h.s. of this equation gives this update 
law properties quite different from a simple gradient search. The functions Xk{t) arising from this 
cost functional may be adjusted to tune convergence properties. In summary, the algorithm works 
as follows: 

1. An initial guess for the control fields u^°\t) is chosen. 

2. The equation of motion (1) with u^°\t) is solved with initial condition g(0) for an ensemble 
of trajectories {g^(t)}. 

3. The co- state ensemble is initialized using the condition A^(T) + M — and then propagated 
backwards to the initial time with Eq. (10). 

4. The ensembles {g%(t)} and {A^(t)} are propagated forward until the end time T, but at each 
time step the co-states are propagated with the old control fields, that is, the ones obtained 
from the j'-th iteration, whereas the states g^(t) are propagated with the new control fields 
given by the update law (11). (Alternatively, stored values from step 3 can be used). 
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5. The objective functional (2) is computed, and steps 3 and 4 are repeated until the value of 
objective functional is in the desired range. 

This optimization algorithm applies not only to Eq. (1), but to any type of linear dynamics 
subject to additive or multiplicative noise. 

III. PROOF OF MONOTONIC CONVERGENCE 

The procedure outlined above guarantees that each step of the iteration reduces the objective 
functional, that is, E[B^ +1 \t), g^ +1) (t)]] < E[B € [u^(t), Q^\t)}}. This becomes evident if the 
change of the objective functional is decomposed into three terms, 



E[L € [uW(t), ^>(t),^}] - E[l 6 [u^(t), g^>(t),^}] = A, + A 2 + A 3 , (12) 



where 



Ai = 



J dtE[R{t,uV\t),$ +1 \t))-R(t,uV\t),$\t))] (13) 

A 2 = ^ T dtE[R(t, u y +1 )(t),^ +1) (t))-R(t,u«(t),4 J+1) (t))] (14) 
G(gf(T))-G($ +1 \T))] (15) 



A 3 = E 

Now Ai > and A3 > follow from Eq. (8) and Eq. (7), respectively, whereas A 2 > is 
a property of the update law (9). We have thus demonstrated monotonicity for the modified 
algorithm, which ensures convergence for an arbitrary initial guess with any objective functional 
bounded from below. 

Recent work by S. Schirmer and P. de Fouquieres [New J. Phys. 13, 073029 (2011)] points 
out that proofs of monotonicity for Krotov-type iteration methods are not rigorous if time-discrete 
approximations are substituted for the dynamics of state and co-state. For all practical purposes, 
however, the typical benefits of the Krotov method persist, i.e., convergence for a wide range of 
initial guesses and rapid improvement during the first iteration steps. 

IV. DRIVEN DYNAMICS WITH FIXED RWA DISSIPATOR: EFFECTIVE RESERVOIR 



For comparison, we study the control problem of the parametrically driven harmonic oscillator, 
modeling the open system dynamics through a Lindblad Master equation. From a system-reservoir 
model, this equation results from performing the Born, Markov, and RWA approximations on the 
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system-reservoir interactions. In the dimensionless units introduced in the main text, the Master 
equation reads 

j f p = -i[H(t),p] 

+ jo(N p + 1) (apa) - \{tfa, p}) 

+ j Np (a! pa - \{aa),p}} (16) 

with 

H(t) = fltfl + \ + ^(a + fit)' = (l + (otfi + I) + ^(a 2 + a*) . (17) 

Following common practice (see, however, Refs. [10, 11]), we have defined the rotating frame 
of the RWA through the time-independent Hamiltonian Hq = a) a + |, not the full Hamiltonian 
H(t). Accordingly, the dissipative terms of Eq. (16) contain the thermal occupation number 
A/3 = (e 13 — l) -1 of the unperturbed oscillator. 

Due to the mismatch between the Hamiltonians of the pre-set rotating frame and of the driven 
dynamics, oscillatory terms reappear when the dissipative terms are transformed to the interaction 
picture. In order to show this, we define the interaction picture with respect to H(t), 

p(t) = U{t)p\{t)U\t) , U = iHU . (18) 

The Master equation for p\ then reads 

j t Pi = io(Np + l)(b0- \{b%pi}) 

+ ^Np^pib-liW,^}) , (19) 

where b — Wall and its adjoint are time dependent, following the equations of motion 

b = — i (! + #)& — *# 6 f (20) 
= if 6 + i(i + f)6t . (21) 

From this, one easily concludes that & is a linear combination of a and fit with [b , w] = 1, i.e., the 
map a I— > b is a Bogoliubov transformation 

b = e ie cosh y a + e^ sinh y fit (22) 

with real (and, in our case, time-dependent) parameters 6, (ft, and y. Using this equation, we obtain 
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the form 

^Pi = 7o(N/3 + 1) cosh 2 y (apia* 1 - i{a f a, pi} ) 

+ 7o(N/3 + 1) sinh 2 y (a* p\a - \{aa\ pi}^j 

+ j Nf3 cosh 2 y (a>pia - \{aa\ pi}) 

+ 7oA> sinh 2 y (ap\a ] - \{a)a, p\}^ 

+ TopNp + l) sinh y cosh ye l(6 ~^ (ap\a - ±{a 2 , Pl }^j 

+ 70(2^ + 1) sinh y cosh ye"^"^ (a t ^ia t - \{a^,pi}) (23) 

for the Master equation, where all of the prefactors depend on time through 9, tp, and y. Typical 
solutions lead to oscillating prefactors, meaning that Eq. (23) can be further simplified through 
time coarse graining. This in itself can be seen as evidence that the initial application of the 
RWA was not consistent. However, coarse-graining at this point is appropriate if one wants to 
compare the features of the dynamics predicted by Eqs. (16) and (23) with the exact dynamics. 
The first four terms of Eq. (23) contain functions of time which are squares of real time-dependent 
functions, therefore they typically dominate the result of the coarse graining procedure. Collecting 
the dominant terms results in the simplified Master equation 

^pi = 7o(^V+ 1) (apia) - \{a)a,pi}^ 

+ joN^ha- liaa^pi}^ (24) 

with 



2N + 1 = (2N p + 1) cosh(2 2/ ) > 2N P + 1 , (25) 

where the bar indicates time coarse graining. The immediate interpretation of this result is obvious: 
Up to a unitary transform, the dynamics of Eqs. (16) and (23) effectively describes relaxation 
towards a thermal oscillator state with a higher-than-thermal occupation number N. This state 
therefore has a higher entropy than the thermal state obtained by equilibrating without driving. 
Accordingly, with the dynamics approximated by Eq. (16), optimal control cannot find solutions 
which push the system entropy below its equilibrium value. This is consistent with our numerical 
findings, in particular the observation that optimal control based on Eq. (16), with overlap to the 
ground state as objective, converges to the trivial result A = 0. The alluringly simple method of 
constructing the dissipator based on a Hamiltonian without driving is an oversimplification which 
fails to reproduce cooperative effects of driving and dissipation, most notably, the cooling effect 
reported in the main text. 



